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I .  INTRODUCTION 


The  Array  has  used  sensitivity  testing  for  many  years,  especially  in  the 
areas  of  primer  testing  and  armor  penetrating  projectiles.  In  sensitivity 
testing,  only  one  observation  per  sample  is  obtained.  This  observation  is 
characterized  by  a  quantal  response,  i.e.,  go  or  no-go,  explosion  or  non¬ 
explosion,  penetration  or  nonpenetration.  The  purpose  of  most  sensitivity 
experiments  is  to  determine  the  amount  of  stimuli  needed  in  order  to  achieve 
a  50%  probability  of  success. 

For  armor  penetrating  projectiles,  the  stimulus  is  velocity.  So  in 
sensitivity  testing,  the  aim  is  t'  determine  the  Vc_,  i.e.,  the  velocity  needed 
to  have  a  50%  probability  of  penetrating  the  armor. 

There  are  many  methods  or  schemes  for  collecting  data  in  sensitivity 
tests.  Some  of  the  more  popular  ones  are  tne  Bruceton  up-and-dovm  method,  the 
Langlie  design,  the  Robbins-Monro  design  and  the  Alexander  design.  For  a 
discussion  of  these  designs  see  Reference  1. 

As  mentioned  earlier,  we  are  interested  in  estimating  and  getting  an 

estimate  of  its  variability.  The  most  commonly  used  technique  of  obtaining 
parameter  estimates  is  to  assume  that  the  critical  stimuli  levels  are  normally 
distributed  and  to  use  maximum  likelihood  estimation  (MLE)  techniques  to 
estimate  y  and  a.  By  critical  stimuli  levels,  i.e.,  critical  penetration 

velocities,  it  is  meant  those  levels  (velocities)  at  or  above  which  a  penetration 
would  always  take  place.  These  estimates  are  not  dependent  upon  the  sensitivity 
design  used  to  collect  the  data. 


There  are  many  computer  programs  that  compute  MLE's. 
these  programs,  at  one  time  or  another,  have  a  converger -j 
initial  values  of  the  estimates  cannot  be  found  that  r 
converging  on  estimates  of  the  mean  critical  stimulus  le 
ard  deviation  of  the  critical  stimuli  levels,  a. 


However,  most  of 
■ .oblem.  That  Is, 
ij  the  prop"  :ji 
,  p,  and  the  star-.i- 


DiDonato  and  Jamagin  used  a  Newton-Raphson  proced«_  ■*  -  ■  .ransformed 

parameters  and  claim  that  it  converges  globally  to  the  "be.  estimates  (when 
estimates  exist),  regardless  of  what  initial  values  were  used.  Unique  MLE's 
will  exist  as  long  as  the  data  meet  the  following  two  restrictions:  1)  min¬ 
imum  level,  a^,  at  which  a  success  is  observed  <  maximum  level,  b ^ ,  at  which 

a  failure  is  observed,  and  2)  average  a^'s  (penetrations)  >  average  b^'s 

''nonpenetrations) .  Tiiey  documented  a  program  written  for  an  IBM  computer 
using  their  procedure.  DiDonato  and  damagin' s  program  was  modified  and 
installed  on  the  Ballistic  Research  Laboratory's  CDC  computer. 


D.  Rothman,  M.J.  Alexander  and  J.M.  Zimmerman,  "The  Design  and  Analysis  of 
Sensitivity  Experiments,"  NASA  CR-62026-Vol.  I,  May  1965. 

^A.R.  DiDonato  and  M.P.  Jamagin,  Jr.,  "Use  of  the  Maximum  Likelihood  Method 
Under  Qviantal  Responses  for  Estimating  the  Parameters  of  a  Normal  Distribution 
and  its  Application  to  an  Armor  Penetration  Prob.’em,"  Naval  Weapons  Laboratory, 
November  1972,  NWL  Technical  Report  TR-2846. 
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This  report  discusses  the  statistical  theory  needed  to  estimate  the  mean 
(Vjq)  and  the  standard  deviation  of  the  response  distribution  and  to  obtain 

the  standard  deviations  of  these  estimates.  A  discussion  of  the  DiDonato  and 
Jarnagin  computing  procedure  used  to  obtain  these  is  provided.  A  nonpararaetric 
algorithm  has  been  incoi^jorated  with  the  DiDonato  and  Jarnagin  program  to 
provide  an  estimate  of  the  mean  (V^q)  and  standard  deviation  of  the  response 

distribucion  when  the  data  do  not  meet  the  requirements  for  the  DiDonato 
and  Jarnagin  procedure.  Examples  are  provided  for  each  of  these  algorithms. 

A  computer  listing,  the  appropriate  inputs  and  corresponding  output  of  tjve  pro¬ 
gram  are  included  in  the  appendices.  The  reader  who  is  more  interested  in  learn¬ 
ing  to  run  the  program  rather  than  in  the  theory  is  referred  directly  to 
Append- X  L. 


II.  THEORY 

principle  of  maximum  likelihood  was  first  introduced  and  extensively 
developed  by  R.  A.  Fisher.  To  aid  in  deriving  the  likelihood  function,  con¬ 
sider  the  results  of  an  experiment  (such  as  an  armor  plate  penetration  problem) 
in  which  the  levels  of  the  stimulus  are  not  completely  under  control.  There  is 
generally  only  one  observation  per  level,  and  the  response  for  each  level  is 
either  a  "success"  or  a  "failure."  Call  the  probability  of  a  successful 
response  p^,  i=l,  ....  n  and  the  probability  of  a  failure  qj ,  j=l,  ...,  m. 

Then  it  follows  from  elementary  statistical  concepts  that  the  probability  of 
the  observed  set  of  responses  may  be  given  by  the  product 

L*  =  PjP2  ... 

Generally,  if  a  total  of  n^m  cases  are  tested  and  there  are  a  successes 
recorded  at  stimulus  levels  a^,  a2,  ...»  and  m  failures  recorded  at  levels 

bj,  b2,  ....  b^,  and  assuming  the  critical  levels  of  stimulus  are  normally  dis¬ 
tributed,  then  the  probability  of  occurrence  of  this  set  of  events  may  be 
written  in  the  more  usual  notation, 

n  m 

L*  =  n  p.  n  q.  (1) 

i=l  ^  j=l  ^ 

(a^-U)/a 

where  p.  =  p[(a.-uVo]  =  I  exp(-t^/2)dt,  (2) 

^  ^  /2?  J 

-00 

and 
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. . . . . . . . 


q,-  =  q[(b. -y)/a]  = 

J  ^  /2? 


exp(-t^/2)dt. 


(bj-y)/a 


The  method  oi  maximuir  likelihood  now  consists  of  choosing,  as  estimates 
of  the  unknown  population  values  of  y  and  o,  those  particular  values  that 
render  L*  as  large  as  possible.  Since  likelihood  functions  are  products,  and 
sums  are  usually  more  convenient  to  work  with  than  products,  the  usual  practice 
is  to  maximize  the  natural  logarithm  of  the  likelihood  function  rather  than  the 
likelihood  function  itself,  i.e.. 


n  m 

Li.  ,0)  =  I  in  p.  +  I  in  q . . 
i=l  ^  j=l  ^ 


The  logarithm  of  the  likelihood  function  has  its  maximum  at  the  same  poin' 
as  does  the  likelihood  function. 

The  maximum- likelihood  estimators  of  y  and  a  are  obtained  by  setting 
the  partial  derivatives  of  the  likelihood  function  equal  to  zero,  i.e.. 


The  likelihood  equations  can  then  be  solved  iteratiyely  using  some  technique 
such  as  Newton-Raphson,  to  obtain  estimates,  y  and  o. 

Previously  it  was  assumed  that  there  was  a  given  set  of  observed  responses, 
n  successes  at  levels  aj^,  a2,  ...,  a^  and  m  failures  at  levels  bj,  b2»  ...,  bj^^. 

Suppose  that  for  a  set  of  distinct  levels  of  stimuli,  there  is  a  statement 
associated  with  each  level  that  it  was  either  a  "success”  or  a  "failure."  Con¬ 
sider  the  discrete  random  variable  (k  =  1,2,...,N),  where  N  equals  the  number 
of  levels  observed  and  each  may  take  on  one  of  two  possible  values,  i.e., 

0,  if  a  failure  occurs  at  level  Z, 

^  = 

I,  if  a  success  occurs  at  level  2.j^. 

If  the  probabilities  of  success  at  the  resulting  levels  are  pj^  and  the 
probabilities  of  failure  are  qj^  =  ^"Pr’  likelihood  of  observing  some 

set  of  responses  (k  =  1,2, ...,N)  at  levels  (k  =  1,2,. ..,N)  is 

L  =  n  p  ^kqJ^'V,  (6) 

k=l  K 


and  the  natural  logarithm  of  the  likelihood  function  becomes, 


(7) 


It  is  this  form  of  the  likelihood  function  that  will  be  used  to  construct 
the  elements  of  the  covariance  matrix  utilizing  the  theory  of  large  sample  dis- 

3  A  A 

tributions.  The  variances  of  y  and  a  may  be  obtained  using  the  relationships 
below  between  the  expected  values  of  the  first  and  second-order  partial  deriva¬ 
tives  of  (Appendix  A), 


I  SL, 


=  E 


3Lo\2 

Sir 


y  ^ 

k=l  Pk‘^k 


N 

I 

k*l 


2  2 
^k  ^k 

Pk^^k 


(8) 


(9) 


and 


a^L 
_ o 

Syao 


at 


aL 


30 


N 


^  k=l  Pk^^k 


(10) 


Additional ly, 


a^L  \ 


-E 


ay' 


‘I 


-E 


ayao 


El— 2.1  =A 

3y 

3L 


yy 


3L. 


"  5irj^53“] 


yo 


(11) 


(12) 


and 


-E 


a^L 


ao' 


31  \2 

^  ‘  j  "  ^0 


(13) 


The  covariance  matrix  for  the  distribution  of  the  maximum  likelihood 

A  A 

estimators  y  and  a  may  be  written  as 


/  V 

-1 

a  1 

\  \a 

^0 

4“° 

A^^j 

(14) 


^A.  Mood  and  F.  Graybill,  Introduction  to  the  Theory  of  Statistics,  McGraw-Hill 
Book  Company,  Inc.,  1963,  pp.  235-239. 
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where  A  and  are  the  asymptotic  variances  of  y  and  o,  respectively,  and 
is  the  asymptotic  covariance  of  y  and  o.  It  should  be  noted  that  a'^^  is 
the  asymptotic  variance  of  the  estimate  of  while  A*^*^  is  the  asymptotic 
variance  of  the  standard  deviation,  a,  of  the  response  distribution. 

For  convenience  of  notation,  equations  (8)  -  (13)  may  be  rewritten. 

Recall  that  k  =  1,2,...,  N,  and  divide  N  into  n  successes  (penetrations)  and 
m  failures  (nonpenetrations) ,  so 

N  =  n  +  m. 

The  n  successes  recorded  at  stimulus  levels  a^,  ^2,  ••.,  a^^  and  m  failures 
recorded  at  levels  bj^,  b2,  ...,  b^  may  be  rewritten  as 

Sj^  =  (a^  -  y)/a,  where  i  =  1,2,...,  n  (15) 

and 

tj  =  (bj  -  y)/a,  where  j  =  1,2,....  ra,  (16) 

which  are  their  standardized  normal  equivalents.  The  s.  are  associated 

with  successes  and  the  t.  are  associated  with  failures.^ 

3 

Also, 

=  z(sp  =  (l/i  (-s^  /2),  when  a  success  is  observed  (17) 

and 

yj  =  z(tj)  =  (l*'^)  exp  •'  ^  /2),  when  a  failure  is  observed.  (18) 

Equations  (8)  -  (10),  which  rej  -jt  the  expected  values  of  the  second-order 
partial  derivatives  taken  ov****  i  entire  set  of  observations,  can  each  be 
partitioned  into  the  sum  ove  i  set  of  successes  plus  the  sum  over  the  set 
of  failures. 

“'2 

The  coefficients  of  a  ,  A  ,  A^^l  may  be  expressed  as 

■jy’  ya 
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o^A  =  y  y 

era  .'s  P-q-  •  t  P'Q* 

1=1  *^1^1  3=1  *^3^3 


III.  PROGRAM  DESCRIPTION 

This  section  first  presents  an  overview  of  the  DiDonato-Jamagin  procedure, 
including:  i)  the  repararaeterization  of  the  likelihood  function,  ii)  the 

A 

actual  procedure  used  to  obtain  p  and  a,  and  iii)  a  description  of  the  output 
sheet.  Should  a  more  detailed  mathematical  description  of  the  analyses  be 
required,  the  reader  is  referred  to  the  DiDonato  and  Jamagin  study. 

A  nonparametric  algorithm  which  has  been  incorporated  into  the  original 
program  and  its  associated  output  sheet  are  also  summarized.  This  procedure 
closely  resembles  a  procedure  presently  utilized  by  the  Materiel  Testing  Direc¬ 
torate,  Aberdeen  Prodng  Ground,  Maiyland. 

The  complete  program  is  available  for  use  on  the  BRL  CDC  computer  system. 
Both  the  DiDonato-Jamagin  procedure  and  the  nonparametric  algorithm  are 
written  in  FORTRAN  IV. 

A.  DiDonato-Jamagin  Procedure 

1.  Likelihood  Function.  Recalling  the  likelihood  function  for  an 
observed  set  of  n  successes  at  stimulus  levels  a.,  a,,  ...»  a  and  m  failures 
at  stimulus  levels  b^,  b2,  ...»  b  ,  we  have 


n  m 

L(P,cr)  =  y  p.  +  y  ^n  q.. 

i=l  ^  3=1  ^ 


DiDonato  and  Jarnagin  then  apply  the  following  transformations  to  replace 
the  parameters  y  and  o  with  new  parameters  a  and  B. 

a  =  y/a  [y  =  a/3]  (23] 

3  =  1/a  >  0  [a  =  1/3].  (24) 

In  addition,  the  original  stimulus  levels  are  reexpressed  as, 

s.  =  a^3-a  (25) 

t.  =  b.3-a.  (26) 

so  that  pj^  and  q^  are  transformed  in  terms  of  the  new  parameters  to 
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The  {a^  =  stimulus  level  at  which  there  is  a  response,  penetration)  and 

{’o.  =  stimulus  level  at  which  there  is  a  nonresponse,  nonpenetration), 

where  i=l,  ...,  n,  j  =  l,  ...,  m,  represent  input  to  the  program.  The  quantities 
a^  and  bj  may  take  on  any  arbitrary  real  values  although  they  are  almost  always 

positive.  The  standard  deviation,  a,  and  its  reciprocal,  3,  are  by  nature 
always  positive. 

In  "heir  study,  DiPonato  and  Jarnagin  point  out  that  to  insure  the 
eXi-tence  of  a  unique  point  at  which  L  attains  a  maximum,  the  a^  and  b.  must 
satisfy  certain  necessary  and  sufficient  conditions,  i.e., 

a„.  <  b  [zone  of  mixed  results]  (37) 

min  max  '■ 

where  a  .  =  min{a.)  and  b  =  max(b.), 

min  ^  1  max  1 


2.  Computing  Procedure.  The  DiDonato- Jarnagin  procedure  first  insures 
that  rhe  necessary  and  sufficient  conditions  indicated  in  (37)  and  (38)  are 
satisfied  by  the  input.  If  either  condition  is  not  satisfied,  then  no 
maximum  likelihood  estimates  exist  for  y  and  a  and  the  program  is  routed  to 
the  nonparametric  algorithm.  If  both  conditions  are  satisfied,  an  initial 
point  (0^,3^)  is  obtained  (unless  the  user  supplies  initial  estimates. 

Appendix  C)  vising  the  following  relationships. 


.  n  m 

where  v  =  -  (72.+  7b.).  (41) 

1=1  3=1 


After  choosing  an  initial  point,  the  ordinary  Newton- Raphson  (N-R)  pro¬ 
cedure  is  applied  in  order  to  reduce  La  and  L3  to  zero  simultaneously  and  produce 
the  N-R  increments  Aa  and  A3  and  a  new  point  (aj,3j).  The  two  iterative  equa¬ 
tions  utilized  by  the  estimation  scheme  are 

AaLaa  +  A3La3  =  -La  (42) 

AaLo3  +  A3L33  =  -L3.  (43) 
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Solving  this  linear  system  of  equations  yields  the  increments  Aa  and  AB 
as  follows: 


Aa  =  (LB  LaB  - 

La  LBB) /A 

(44) 

AB  =  (La  LaB  - 

LB  Loa) /A 

(45) 

where 

A  =  Loa  LBB  - 

L^aB  >  0. 

(46) 

The  increments  Aa  and  AB  are  the  changes  in  the  old  values  of  a  ard  B  calcu¬ 
lated  at  each  iteration  of  the  N-R  procedure. 

If  a  sequence,  (a.  ,B^),  determined  by  iterativjn  and  the 

rh  ^  /V  ^ 

and  Bj^  denote  the  approxi.nations  to  a  and  B,  then  the  iterations  are 
terminated  when 

lAOkl  <  and  jABj^l  < 

The  parameters  and  ^2  defined  as  follows: 

Ej  =  2.5  X  i0~^  =  1/2  £.2' 

At  this  point,  additional  notation  is  introduced  in  order  to  rewrite  'equations 
(19)  -  (21)  and  (30)  -  (34)  as  actually  used  in  the  DiDonato-Jamagin  procedure. 

Referring  to  the  parameter  input  description  in  Appendix  C,  it  can  be 
seen  that  only  those  a^^  which  are  different  are  listed  along  with  an  integer 

n(i)  which  denotes  the  number  of  times  a^  appears  as  input.  The  integer 

m(j)  denotes  the  number  of  times  each  different  b^  appears  as  input.  Noting 

that  the  expressions  for  La,  LB,  Laa,  LoB,  LBB  are  linear  sums,  quantities 

such  as  (Xj^/pj^)  and  (y^/q^)  need  only  be  computed  for  the  different  a.  or  b^ , 

then  multiplied  by  n(i)  or  m(j),  respectively. 

If  the  k^^  different  a.  is  denoted  by  a(k),  the  r^^  different  b^  is 
denoted  by  b(r),  and  K  and  R  denote  the  total  number  of  different  a^and  b^, 
respectively,  then 

K  R 

n  =  f  n(k),  m  =  f  m(r)  (49) 

k=l  r*l 

and  n  and  ra  have  their  usual  meaning.  The  equations  referenced  above  can  now 
be  written  as: 


(47, 

(48) 
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(50) 


La  =  , 
r: 


LB  = 

k 


Low  = 


Lap  = 


LBB  = 


"'2 

o  Ayy  = 


‘'2 

o  Amo 


[  in(r)Cy-/0  -  I 

=1  k=l 

K  ^ 

I  n(k)a(k)(x^/Pj^)  -  ^I^m(r)b(r)(y^/qp 

R 

I  raiT)iy^/q^)(.y^/%  - 

T=1 

I  n(l‘)(\/Pk)K/Pk  ^ 

c=l 

K 

y  n(k)a(k)CA,^/p^)(x^/Pj^  +  s^) 

k=l 

R 

y  m(r)bCr)  (y^/qp  (y^Qj.  "  ^P 
r=l 

R  ■> 

-  y  ra(r)b  (r)  (ypqp  (yp^y  "  ^P 
r=l 

K  , 

-  y  n(k)a^(k)(xpp,^)(xpp^  +  sp 

k=l 

f  n(k)(xppp(xpqp 
k=l 

R 

y  ni(r)(yppp(ypqp 

r=l 

K 

y  n(k)s^(x^/pp(x^/qk) 

k=l 

R 

^  y  in(r)t^(yppp(ypqp 

r=l 


(51) 


(52) 


(53) 


(54) 


(55) 


(56) 
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'  alphanumeric  characters  which  may  indicate  a  title,  identification  number,  etc. 

This  is  followed  by  a  listing  of  the  different  a.  and  b.  and  the  number  of 

X  3  A.  /w 

times  each  value  occurs.  The  maximum  likelihood  estimates,  y  and  a,  identi- 
‘  fied  as  'MU'  and  'SIGMA,*  respectively,  are  listed  toward  the  top  center. 

'MU,'  or  y,  is  the  estimate  of  while  'SIGMA,'  or  a,  is  the  estimate  of 

1  the  standard  deviation  of  the  response  distribution.  Directly  below  the 

estimates,  the  elements  of  the  covariance  matrix  are  recorded,  followed  by 

A  A 

and  identified  as  the  'STD  DEV  MU*  and  'STD  DEV  SIGMA,*  respectively. 

i  These  quantities  are  followed  by  the  initial  approximations  of  ('ALPH0'), 

Bg  ('BETA0'),  y^  ('MU0'),  and  (’SIGMA0*).  The  Newton -Raphson  increments 

(  at  each  iteration,  Aa,  AB,  written  as  'DELTA  ALPHA'  and  'DELTA  BETA,*  as 

well  as  the  associated  value  of  L*  are  listed  next.  The  last  line  lists  the 
ff.  •’.1  estimates  of  the  likelihood  function  L*,  the  determinant  A,  a,  and  B 
ano  are  identified  as  'MAXIMUM,'  'DELTA,'  'ALPHA,*  and  'BETA,'  respectively. 


B.  Nonparametric  Algorithm 

1.  Computing  Algorithm.  As  previously  mentioned,  this  algorithm  is 
initiated  should  either  of  the  necessary  and  sufficient  conditions  imposed 
upon  the  input  data  not  be  satisfied.  This  routine  will  compute  a  maximum 
of  three  (3)  estimates  for  the  mean  and  standard  deviation  using  a  method 
of  averaging  particular  levels  of  response  and  nonresponse.  To  facilitate 
execution  of  the  algorithm,  the  original  input  arrays  have  been  expanded  to 
include  all  a^'s  and  b^'s,  including  replicated  values.  Each  array  is  rank 

ordered  from  the  smallest  to  the  largest  level. 

The  first  estimate  of  the  mean  is  obtained  by  averaging  the  lowest  value 
of  a  response  and  the  highest  value  of  a  nonresponse.  Then  the  two  lowest 
responses  and  the  two  jargesi  nonresponses  are  averaged  to  obtain  a  second 
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estimate  for  the  mean.  This  iterative  procedure  may  be  continued  until  all 
even  number  combinations  of  the  a^^'s  and  b^'s  are  exhausted.  However,  it 

has  been  elected  to  allow  the  algorithm  to  generate  only  a  possible  maximum 
of  three  (3)  estimates  (requiring  a  total  of  six  (6)  input  levels).  At  each 
iteration  an  estimate  of  the  standard  deviation  of  the  distribution  is 


calculated  as  follows: 

1/2  ^  [range  between  levels]  (60) 

where  'range  between  levels'  is  defined  as  the  difference  between 

max  [B0(MTOTB).  A0(I)J  -  rain  [B0(J),  A0(1)]  (61) 

and  A0(1)  =  smallest  response  value  of  input  array  (62) 

A0(I)  =  response  value  at  iteration,  k=l . 3  (63) 


B0(J)  =  nonresponse  value  at  k^''  iteration,  k=l,  ...,3(64) 


B0(MTOTB)  *  largest  nonresponse  value  of  input  array.  (65) 

2.  Output  Description.  The  output  generated  for  example  2  (Section  4) 

is  shown  in  Appendix  D.  At  the  top  of  the  output  sheet,  the  program  prints 
out  a  description  of  the  first  necessary  and  sufficient  condition  that  is 
not  satisfied.  Directly  under  this  is  an  ordered  listing  of  the  input  arrays 
of  responses  and  nonresponses,  including  replications.  Finally,  for  each 
iteration,  the  estimate  number  and  the  number  of  rounds  required  for  each 

A,  A 

estimate  of  y  and  o  followed  by  the  values  for  y  and  o,  written  as  'MU'  and 
'SIGMA,'  respectively,  are  shown. 


IV.  APPLICATIONS  (EXAMPLES) 

The  three  cases  used  for  illustration  ai*e  taken  from  armor  plate  pene¬ 
tration  studies.  The  actual  output  of  the  computer  program  for  each  of  the 
examples  is  included  as  Appendix  D.  The  data  set  for  the  first  example 
satisfied  the  necessary  and  sufficient  conditions  outlined  in  Section  III; 
therefore,  the  output  will  be  generated  from  the  DiDonato-J.imagin  procedure. 
Example  2  is  a  'no  zone  of  mixed  results'  case  which  utilized  the  nonparametric 
algorithm.  The  third  example  satisfied  both  conditions  required  for  MLE's 

A 

to  exist  for  y  and  o,  yet  diverged  using  computer  generated  starting  values 
in  a  different  estimation  program.  This  case  is  shown  for  two  different 
sets  of  starting  values  for  which  the  results  shew  convergence. 
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For  each  example  the  initial  estimates  were  generated  by  the  program; 
in  addition,  the  third  example  was  run  a  second  time  with  user  supplied  start¬ 
ing  values  for  a  and  B  . 

o  o 

EXAMPLE  1 

in  firing  10  rounds  at  a  given  armor  plate  the  following  observations 
were  recorded: 


VELOCITY  (ro/s) 


RESPONSE 

Penetration 

Nonpenetration 

Penetration 

Penetration 

Penetration 

Penetration 

Nonpenetration 

Nonpenetration 

Penetration 

Nonpenetration 


Initial  estimates  computed  were  =  37.567  m/s  =  956.625  m/s), 

B  -  .039  m/s  (o^  =  25.464  m/s)  The  N-R  iteration  procedure  ^as  then  applied 

4  times  to  generate  final  estimates  of  a  =  32.123  m/s  (p  948.823  m/s), 

A  A 

B  =  .034  m/s  (o  s  29.537  m/s).  Employing  the  solutions  obtained,  the 

A  A 

approximate  asymptotic,  variances  of  p  and  o  were  obtained.  The  covariance 
matrix  is  given  by 


183.506 


-  50.653 


-50.653 


358.946  . 


Thus,o.‘'^  =  183.506  ra/s  {a"  =  13.546  m/s)  and  =  358.946  m/s  (o"  =  18.946  m/s), 
y  u  0  0 

EXAMPLE  2 

Seven  rounds  of  a  given  projectile  were  fired  at  a  giver,  armor  plate.  The 
following  observations  were  recorded: 


tiliwillilfeM 


VELOCITY  (m/s) 

RESPONSE 

944 

Nonpenetration 

973 

Penetration 

961 

Nonpi  \etration 

982 

Penetration 

966 

Penetration 

949 

Nonpenetration 

970 

Penetration 

Once  the  data  were  sorted,  it  was  clearly  seen  that  there  was  'no  zone  of  mixed 
results.'  The  program  then  routed  to  the  nonparametric  algorithm.  To  illus- 

A  A 

trate  the  method,  I^t's  calculate  the  second  estimate  of  HU  (p)  and  SIGMA  (o) . 
The  second  estimate  required  a  total  of  4  levels  or  rounds  {2  lowest  responses, 
2  highest  nonresponses).  From  the  sorted  input  lists,  the  two  lowest  responses 
were  966  m/s  and  970  m/s;  the  two  highest  nonresponses  were  949  m/s  and 
961  m/s.  From  the  routine, 

WJ0  =  MU0  ♦  (A(2)  ♦  B(2j)/2.0 

Lstimate  1  for  MU  =  963.5,  so 

MU0  =  963.5  ♦  (970  +  949)/2.0 

MU0  =  1923 


W  =  2.0  (W  =  KOUNT;  KOUNT  =  iteration  ...• 
MU  =  MU0/W  =  1923/2.0  =  961.5  m/s 
SIGMA  =  (RRANGE  -  SRANGE)/2.0 
Where  RRANGE  =  MAX(B0(MTOTl"»  ,  A0(I)) 

=  MAX (961,  970) 

=  970 

SRANGE  =  MIN(B0(J),  A0(1)) 

=  MIN (949,  966) 

=  949 

Thus,  SIGMA  =  (970  -  949) /2.0 
=  10.5  m/s 
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EXAMPLE  3 


Again,  10  rounds  were  fired  at  a  given  amor  plate  and  yielded  the 
following  results: 


VELOCITY  (m/s) 

RESPCMSE 

1295 

Nonpenetration 

1296 

Penetration 

1298 

Nonpenetration 

1301 

Nonpenetrat ion 

1303 

Nonpenetration 

1304 

Nonpenetration 

1307 

Nonpenetration 

1310 

Penetration 

1311 

Penetration 

1314 

Nonpenetration 

a)  Initial  estimates  generated  by  the  program  were  »  210.224  m/s 
(p  =  1304.405  m/s),  =  .161  m/s  {o  =  6.205  m/s).  Applying  the  N-R  scheme 

0  O  aO  a 

3  times  produced  final  estimates  of  o  =  50.657  m/s  (y  *  1317.893  m/s), 

A  A 

3  *  .038  m/s  (0  =  26.016  ro/s).  The  following  asymptotic  variances  were 
obtained : 

oj  =  690.476  m/s  (oj)  =  26,277  m/s) 
and 

a|  =  2142.252  ra/s  (o^  =  46.283  m/s). 


b)  Initial  approximations  of  *  26.000  m/s,  =  .006  m/s  were  input 
by  the  user.  These  values  generated  initial  values  for  and  of 
4333.333  m/s  and  166.667  m/s,  respectively.  The  N-R  procedure  had  to  be  applied 

A  A  A  ^ 

5  times.  Final  estimates  for  o,  B,  p,  o  were  50.658  a/s,  .038  m/s, 

1317.892  m/s,  and  26.016  m/s,  respectively.  Applying  these  solutions  yielded 

A  A 

the  following  asymptotic  variances  for  p  and  o; 


V.  SUNWARY 


For  problems  of  a  quantal  response  nature,  a  discussion  of  the  maximum 
likelihood  estimation  theory  and  large-sample  distribution  theory  required 
to  obtain  the  estimates  of  the  mean  (Vr^)  and  standard  dc/iation  of  a  cumula- 

DU 

tive  normal  response  function  has  been  g'ven. 

The  A.  R.  DiUonato  and  M.  P.  Jamagin,  Jr.  maximum  likelihood  estimation 
program  has  been  modified  and  installed  on  the  Ballistic  Research  Laboratory's 
CDC  computer.  Instructions  needed  to  properly  execute  the  program  have  been 
provided.  This  program  has  successfully  converged  to  provide  maximum  likeli¬ 
hood  estimates  for  all  data  sets  tried  by  the  authors.  A  nonparametric 
algorithm  is  also  available  for  those  cases  where  maximum  likelihood  estimates 
do  not  exist. 
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APPENDIX  A 


EXPECTED  VALDES  NEEDED  FOR  ASYMPTOTIC  VARIANCES  OF  p  AND  0 
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. . am 


3L  N 


sr  °  ‘  \  Sir  Pk  *  3u  \  > 


N  ,  -z.  ,  z, 

=  I  [  •  —  •  -7;^  +  •  i-  ‘  ] 

k=l  ^  h  ^  ° 


,  N 
0  1  r 


ST  =  j  t  »-^k) 


N  - 

o  _  r  r  r  9 


sr  “  [  ''k  3j  Pk  *  »-*k'i  to  Pk  i 


=  J  1  6^  •  i  •  .  (1.5J  .  j-  .  ) 

kil  k  Pk  «  ‘  k'  0  ' 


1  H  ^th  ''ich 


Now,  consider  the  fact  that  ‘^[g  =  T  gCS^j'  ,,)  ,  where 

k=l 

E(6j^)  -  f(l)  =  pj^  and  E(l-6|^)  =  f(0)  =  qj^;  then  the  expected  values  of  equa¬ 
tions  (A5)  and  (A6]  become 


3L  ,  N 
0.1  r 


~  n^K  ■ 

o  k=I 


3^L 


-E  [ 


2  „  2 
,  N  v.z.  N  v,z. 

]  =  Vl  I  4^  vi 


3y3a  ^2 


k=l  Pk 


Therefore,  equations  (A9)  -  (All)  may  be  rewritten  as 


and 


a^L 


U  LI  N  z* 

-E  t  — I  1  •  ^  I  ^  ' 

3p  a  k=l  »'k^k 


3^L. 


-E  [ 


3a 


M  2  2 

1  =  V  I 


2  ^  ■  a2  p,q,  ' 


3^L 


-E  f  ^1  =-2  .1 


^  ''k^k 


a^  k=l  Pk'lk 
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(All) 


(A12) 


(A13) 


(A14) 
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APPENDIX  B.  PARTIAL  DERIVATIVES  NEEDED  FOR  NEWT(»i-RAPHSON  PROCEDURE 


We  wish  to  compute  the  partial  derivatives  of  the  p^  and  as  an  aid 
in  computing  the  partial  derivatives  of  the  function. 


n  m 


and 


3^  =  exp[-(a.e-a)2/2]  =  a^x. . 

2  2 

3  p.  -a.  (a.g-o)  2 

^  ^  ^  -  exp[.(a.B-a)72l 


36 


3  Pi  2 


_  =  -a.  s.x. . 
2  111 


36 


Similarly  for  the  q^. 


Let  3q.  2 

=  (1/^)  exp[-(bj6-a)V2]  = 


then,  3  q.  3y.  (b.6-a)  2 

- i  =  — i-  =  — i- -  exp[-(b.6-o)  /2] 

3a^  3a  ^ 


3^q, 


3a 


t  =  t,y 


rr 


3  q.  3y.  -b.(b.6-a)  2 

- 1  =  — i  =  — I—il - exp[-(b.6-a)  /2] 

3a36  36  ^ 


3o36  ' 


^  exp[-(b.6-a)^/2]  =  -  b.y., 
36  ^  J  1  J 


(B9) 


(BIO) 


(BID 


(B12) 


(B13) 


(B14) 


and  3^q.  b.^(b.6-a)  2 

— ^  =  J  -  .-j- - exp[-(b.8-o)  /21 

3g2  3 


2 

=  b,  t.y,. 
36 


.2  j  ri 
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(BIS) 


Using  che  above  results,  we  have 


La 


'  H  =  .1 5I fO  •  j,  fe  "li 


i=l 


?  J.  !!i.  ? 

■  i';i  n  3=1  ’) 

La  •  I  ■  J, 

j=l  J  •'  i-1 


=  I 

3=1 


3y.  SQj 


n  Pit  TS±2j_!iLlili 

I  2  “ 

i=l  Pi 


n>  t.y.  y 


n  -s.x.  X 

2  ' 


»u  '’A/i  ^  r  /  A 

3=1  ^3  ‘Ij  ^ 


Loa 


.  -  !  (yi/q,)(Cyj/q,)-Ljl  -  j, 

>1  3  J  3  J  J  1=1 


(Bib) 


,  (B17) 
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APPENDIX  C.  INPUT  NEEDED  AND  PROGRAM  LISTING 

The  A.  R.  DiDonato  and  M.  P.  Jamagin,  Jr.  maximum  likelihood  estimating 
program  is  available  on  the  Ballistic  Research  Laboratory's  CDC  computer 
system. 

Prior  to  executing  the  program,  a  data  file  must  be  created  and  stored 
on  MFA  (mainframe  A) .  A  list  of  the  parameters  that  must  be  read  in  as  input 
and  the  format  and  instructions  needed  to  create  the  data  file  are  provided 
here. 


The  three  (3)  card  jobstream  that  must  be  created  and  submitted  to  MFZ 
(mainframe  Z)  in  order  to  run  the  program  is  also  given. 

It  is  important  to  remember  that  the  data  file  must  be  made  accessible 
to  the  MFZ  while  in  the  interactive  facility  (lAF)  mode. 

The  following  list  briefly  describes  each  of  the  parameters  which  must 
be  read  into  the  program  prior  to  execution  and  which  also  comprises  the 
calling  sequence  for  the  main  computing  subroutine  EPPA. 


CALL  EPPA  (IDENT,K,L,ALPHO,BETO,FNA,A,FNB,B) 

IDENT  -  Array  dimensioned  at  8  locations  with  up  to  80  characters 
allowed  per  job.  The  Hollerith  character  content  of  IDENT 
is  printed  on  the  top  line  of  the  output.  Note:  IDENT  is  not 
printed  when  the  nonparametric  procedure  is  employed. 

K-  Twice  n,  where  n  is  the  number  of  numerically  different  a. 's 
(responses) .  n  >  2. 

L-  Twice  m,  where  m  is  the  number  of  numerically  different  b.’s 
(nonresponses).  m  >  2.  ^ 

ALPHO.BETO  -  User  supplied  starting  values  o  ,  respectively,  to  the 

algorithm.  Setting  BETO  <  0  allows  ?he  user  to  have  the  algorithm 
compute  the  starting  values  ®o  ” 

FNA, A  -  Arrays  dimensioned  at  k  =  2»n.  FNA(i)  specifies  the  number  of 

A(i)  values  used,  i  =  1,2,  ...,n.  FNA(n+l),  FNA(n+2),  ..., 

FNA(2n)  and  A(n^l),  A(n'*-2) ,  ...,  A(2n)  is  used  by  EPPA  as 
workspace. 

FNB, B  -  Arrays  dimensioned  at  £  =  2»m.  FNB(i)  specifies  the  number  of 

B(i)  values  used,  i  =  1,2,  . . .  ,m.  ^B(m'»-1),  FNB(m'*-2),  ...FNB(2m) 

and  B(m+1) ,  B(m+2),  ...,  B(2m)  is  used  by  EPPA  as  workspace. 
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STEP  1.  CREATING  THE  DATA  FILE 


Data  should  be  created  and  stored  as  an  MFA  permanent  file. 


Card  1 

Columns 

Format 

IDENT 

1-80 

8A10 

Card  2  2  inputs 

K 

1-5 

15 

L 

6-10 

15 

Card  3  2  inputs 

ALPHO 

1-10 

FIO. 2 

BETO 

11-20 

FIO. 2 

Cards  4,5 

Input  may  consist  of  a  maximum  of  8  data  points  per  card.  More 
than  one  data  card  may  be  required  to  input  each  parameter.  These 
cards  contain  data  read  in  with  an  F10.2  format. 


Columns 

1-10 

11-20  . 

Card  4 

FNA(l) 

FNA(2) 

FNA(8) 

Card  5 

Ad) 

A(2) 

A{8) 

Cards  6,7 

Input  may  consist  of  a  maximum  of  8  data  points  per  card.  More 
than  one  data  card  may  be  required  to  input  each  parameter.  These 
cards  contain  data  read  in  with  an  FIO.?  format. 


Columns 

1-10 

11-20  . 

Card  6 

FNBfl) 

FNB(2) 

FNB(8) 

Card  7 

B(i) 

B(2) 

B(8) 

STEP  2. 

Prior  to  running  the  program,  the  data  file  must  be  made  accessible  to 
the  MFZ  while  in  lAF  (/)  mode  with 

PERMIT,  PFN,  MFZ 

where 

PFN  =  filename  under  which  the  input  data  is  created 


4*' 


STEP  3. 


To  run  the  programj  create  and  submit  the  following  3  card  jobstrean  to  MFZ. 

JOBN/tME,  STMFZ. 

ACCOUNT,  XXXXXX. 

BEGIN,  NMLEX,  WiCEX,  PFI= _ ,  UN= _ ,  (RJE=RJE _ ). 

where 

PFI  =  file  name  under  which  the  input  data  is  stored, 

UN  =  user  name  identification  for  above  file, 

RJE  =  RJEXXXX  where  XXXX  is  a  4- digit  code  designating  a  particular 
RJE  as  the  output  device;  if  omitted,  the  central  site  serves 
as  the  output  destination. 
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OUTPUT  FOR  EXAMPLE  3b 


DISTRIBUTION  LIST 


No.  of  No.  of 

Copies  Organization  Copies  Organization 

12  Administrator  1  Commander 

Defense  Technical  Info  Center  US  Aimiy  Communications  Rsch 

ATTN:  DTIC-DDA  and  Development  command 

Cameron  Station  ATTN:  DRDCO-PPA-SA 

Alexandria,  VA  22314  Fort  Monmouth,  NJ  07703 


1  Commander  1  Commander 

US  Army  Materiel  Development  US  Army  Electronics  Research 

and  Readiness  Command  and  Development  Command 

ATTN:  DRCDMD-ST  Technical  Support  Activity 

5001  Eisenhower  Avenue  ATTN:  DELSD-L 

Alexandria,  VA  22333  Fort  Monmouth,  NJ  07703 


1  Commander 

US  Army  Armament  Research 
and  Development  Command 
ATTN :  DRDAR-TDC 
Dover,  NJ  07801 

2  Commander 

US  Army  Armament  Research 
and  Development  Command 
ATTN:  DRDAR-TSS 
Dover,  NJ  07801 

1  Commander 

US  Army  Armament  Materiel 
Readiness  Command 
ATTN:  DRSAR-LEP-L 
Rock  Island,  IL  61299 

1  Director 

US  Army  ARRADCOM 
Benet  Weapons  Laboratory 
ATTN:  DRDAR-LCB-TL 
Watervliet,  NY  12189 

1  Commander 

US  Army  Aviation  Research 
and  Development  Command 
ATTN:  DRDAV-F. 

4300  Goodfellow  Blvd. 

St.  Louis,  MO  63120 

1  Director 

US  Army  Air  Mobility  Research 
and  Development  Laboratory 
Ames  Researcii  Center 
Moffett  Field,  CA  94035 


1  Commander 

US  Army  Missile  Command 

ATTN:  DRSMI-R 

Redstone  Arsenal,  AL  35898 

1  Commander 

US  Army  Missile  Command 
ATTN:  DRSMI-YDL 
Redstone  Arsenal,  AL  35898 

1  Commander 

US  Army  Tank  Automotive  Rsch 
and  Development  Command 
ATTN :  DRDTA-UL 
Warren,  MI  48090 

1  Director 

US  Army  TRADOC  Systems 
Analysis  Activity 
ATTN:  ATAA-SL 
White  Sands  Missile  Range 
NM  88002 

1  Commander 

US  Army  Research  Office 

ATTN:  DRXRO-MA,  Dr.  Robert  Launer 

P.O.  Box  12211 

Research  Triangle  Park,  NC  27709 
1  Commander 

US  Army  Combat  Development 
Experimentation  Command 
ATTN;  ATEC-SA,  Dr.  Marion  R.  Bryson 
Fort  Ord,  CA  93941 
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distribution  list 
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Organization 


2  Commandant 

US  Army  Infantry  School 
ATTN;  ATSH-CD-CSO-OR 
Fort  Benninft,  GA  31905 


Aberdeen  Proving  ~.ouiid 

Dir.  USAMSAA 

ATTN:  DRXSY-D 

DRXSY-MP,  H.  Cohen 

Cdr,  USATECOM 

ATTN:  DRSTE-TO-F 
Dir,  USACSL,  Bldg.  E3516,  EA 
ATTN:  DRDAR-CLB-PA 
DRDAR-CLN 
DRDAR-CLJ-L 


USER  EVALUATION  OF  REPORT 
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